load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/substring.ncl"
load "/public/users/wrf/Codes/WRF_scripts/utils/WRF_dyy_contributed.ncl"


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;            User Config                   ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; radius to degree
r2d = 45.0/atan(1.0)

if (.not. isvar("configFile")) then
    configFile = "/public/users/wrf/Data/NorthEast_gridwind/grid_list_d01.csv"
    print("configFile is undefined, using default: " + configFile)
else
    print("Using configFile: " + configFile)
end if

grid_indices = readAsciiTable(configFile,2,"integer",0)
dim_indices  = dimsizes(grid_indices)
ngrids       = dim_indices(0)

levels             = (/10, 30, 50, 70/)
;levels             = (/10, 30, 50, 70, 80, 120/)
nlevels    = dimsizes(levels)

;;;;;;;;;;;;;;;;;;;;;;
; Env and Settings   ;
;;;;;;;;;;;;;;;;;;;;;;


nlevel = dimsizes(levels)

if (.not. isvar("output_prefix")) then
    output_prefix = "~/wrfpost/test/select_gridwind"
    print("output_prefix is undefined, using default = " + output_prefix)
else
    print("Using output_prefix: " + output_prefix)
end if

if (.not. isfilepresent(output_prefix)) then
    system("mkdir " + output_prefix)
end if

; test if input_file is given
if (.not. isvar("input_file")) then
    input_file = "/public/users/wrf/wrfout/2012120608/wrfout_d01_2012-12-06_00:00:00"
    print("input_file is undefined, using default = " + input_file)
else
    print("Using input_file: : " + input_file)
end if

; default output time interval is 15min
if (.not. isvar("output_interval")) then
    output_interval = 0.25
    print("output_interval is undefined, using default = " + output_interval)
else
    print("Using output_interval = " + output_interval)
end if

a = addfile(input_file + ".nc", "r")
filename = systemfunc("basename " + input_file)

; check start_date
if (.not. isvar("start_date")) then
    start_date = substring(filename, 11, 23)
    print("start_date is undefined, using default = " + start_date)
else
    print("Using start_date: " + start_date)
end if

domain = substring(filename, 7, 9)

times = a->Times
hours = WRF_Times2Udunits_c(times, 0)
ntimes = dimsizes(hours)

; UTC -> CST
hours = hours + 8 ; 
timeStr_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")

wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)
dt = round(output_interval / wrf_dt, 3)

; create time series in desired time resolution
output_hours = fspan(hours(0), hours(ntimes-1), tointeger((hours(ntimes-1) - hours(0))/output_interval + 1))
output_ntime = dimsizes(output_hours)
copy_VarAtts(hours, output_hours)

;;;;;;;;;;;;;;;;;;;;;;;;;
; Loop for each station ;
;;;;;;;;;;;;;;;;;;;;;;;;;

do igrid = 0, ngrids - 1
    i = grid_indices(igrid, 0)
    j = grid_indices(igrid, 1)
    print("Processing grid:(" + i + "," + j + ")")
    if (.not. isStrSubset(configFile, domain)) then
        print("The domain is not desired, skipping.")
        continue
    end if


    output_filename = output_prefix + "/gridwind_" + "gn_" + sprinti("%03i",i) + sprinti("%03i",j) +  "_" + domain +  "_" + start_date + ".csv"


    line_array = new(ntimes(0)*nlevel, "string")
    i_line = 0

    ;;; Load Terrain Data ;;;
    yy = (/j,j+1/)
    xx = (/i,i+1/)
    kk = (/0, 1, 2/)

    ;;; load variables
    ; read variables from wrfout ncfile
    tmp_ter         = wrf_dyy_getvar(a,"ter",0,kk,yy,xx)
    tmp_z           = wrf_dyy_getvar(a,"z",-1,kk,yy,xx)
    tmp_pres_hpa    = wrf_dyy_getvar(a, "pressure", -1, kk, yy, xx)
    tmp_uvmet = wrf_dyy_getvar(a, "uvmet", -1, kk, yy, xx)
    tmp_uvmet10 = wrf_dyy_getvar(a, "uvmet10", -1, kk, yy, xx)

    ter = tmp_ter(0,0)
    z   = tmp_z(:,:,0,0)
    pres_hpa    = tmp_pres_hpa(:,:,0,0)
    uvmet       = tmp_uvmet(:,:,:,0,0)
    uvmet10     = tmp_uvmet10(:,:,0,0)

    ; interpolate to the site location

    ; get surface pressure
    surf_pres = pres_hpa(:,0)

    ; do the vertical interpolation
    uv_arr = wrf_uv_intrp_vert(ter, z, uvmet, uvmet10, levels)

    ; calculate wind speed and direction
    spd_transpose = transpose(sqrt(uv_arr(0,:,:)^2 + uv_arr(1,:,:)^2))
    dir_transpose = transpose(round(atan2(uv_arr(0,:,:), uv_arr(1,:,:)) * r2d + 180, 3))

    ; interpolate to 15 min
    output_spd = dyy_interp_ts(hours, spd_transpose, output_hours, "cubic")
    output_dir = tointeger(dyy_interp_ts(hours, spd_transpose, output_hours, "direction"))
    delete(spd_transpose)
    delete(dir_transpose)
    delete(uv_arr)

    ;;; END OF WIND EXTRACTION ;;;


    ;;; write to txt file ;;;
    header = "I, J,   TIME"
    ; station number and time
    printArray = i + "," + j + "," + ut_string(output_hours, "%Y-%N-%D_%H:%M:%S")
    
    ; wind speed and direction
    do ilvl = 0, nlevel - 1
        header = header + ", " + sprinti("%3i", levels(ilvl)) + "sp, " + sprinti("%3i", levels(ilvl)) + "dir"
        printArray = printArray + ", " + sprintf("%5.2f",output_spd(ilvl,:)) + ", " + sprinti("%3i", output_dir(ilvl,:))
    end do


    allPrintArray = array_append_record(header, printArray, 0)

    asciiwrite(output_filename, allPrintArray)
    print("Creating: " + output_filename)

    delete(ter)
    delete(z)
    delete(pres_hpa)
    delete(printArray)
    delete(allPrintArray)
end do
